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ABSTRACT 

A procedure is developed to determine approximate periodic solutions of autonomous and 
non-autonomous systems. The trignometic collocation method (TCM) is formalized to allow for 
the analysis analysis of relatively small order systems directly in physical coordinates. The TCM is 
extended to large order systems by utilizing modal analysis in a component mode synthesis 
strategy. The procedure was coded and verified by several check cases. Numerical results for two 
small order mechanical systems and one large order rotor dynamic system are presented. The 
method allows for the possibility of approximating periodic responses for large order forced and 
self excited nonlinear systems. 


NOMENCLATURE 


a o 

Fourier static coef. vector 

a i’bj 

Fourier cosine, sine coef. vector 

A,£ 

state matrices 

c 

vector of Fourier coef. 

£ 

matrix of collocation values 

D 

tridiagonal matrix 

f 

state function vector 

F 

state force vector 

£ 

damping matrix 

K 

stiffness matrix 

M 

mass matrix 

m 

no. of harmonics 

N 

no. of collocation points 

n 

no. of modes 

A 

n 

truncated no. of modes 

q 

physical coordinate vector 

q b 

nonlinear subvector of q 

Q s 

linear system force vector 

Q b 

nonlinear force vector 

R 

norm 

r,I 

right, left displacement vectors 

I 

right modal matrix 

£ 

connectivity matrix 

t 

time 

T 

fundamental period 
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transformation matrix 

state vector 

state right, left vectors 


I 

x 

y.z 

GREEK 

8jj Kronecker delta 

X eigenvalue 

r\ modal coordinate 

co fundamental frequency 

SUPERSCRIPT 

o collocation values 

• d/dt 

T transpose 


INTRODUCTION 

Nonlinear phenomena of many forms are clearly present in all complicated machinery. The 
future development and advancements in these machines depends strongly on our ability to 
identify, understand, model, analyze and design with these various nonlinear mechanisms present. 
Transient and steady state analysis capabilities are required with direct numerical integration, 
presently the most popular tool. This work presents a method based on trignometric collocation for 
approximating periodic solutions of forced systems and for locating limit cycles of self excited 
systems. The use of modal analysis allows the method to be extended to large order systems. 

Previous work on the steady state response of systems which include nonlinear 
components is limited except by direct numerical integration. This can be very time consuming, 
especially for large order systems, and is not particularly economical in parametric design 
applications. It is really the only option available for transient analysis, however, and also serves 
as a useful means for verifying final designs. 

Some quantitative methods for steady state analysis of nonlinear systems include 
perturbation techniques, describing function procedures, harmonic balance procedures, and 
methods of weighted residuals. Perturbation techniques (Nayfeh, 1981) have a limited range of 
applicability due primarily to high algebraic complexity for large order systems. They also require 
the introduction of a small parameter, thus restricting the solution validity to systems with weak 
nonlinearities. Describing function methods (Atherton, 1982) are a good choice for many problems 
since they can accommodate non-analytic nonlinearities. They can be used, however, only when 
higher harmonics are small compared to the fundamental component. 

The harmonic balance method, (Hagedom, 1982), has been recently applied to the analysis 
of engineering systems ( Yamauchi, 1983; Saito, 1985) and the preliminary results indicate that the 
method may be quite effective. An alternate approach is the use of methods of weighted residuals 
which have been used quite extensively in the past to solve nonlinear boundary value problems. 
Some of these methods, which have been extended to the problem of determining periodic 
response, include Galerkin’s method (Urabe, 1965; Urabe and Reiter, 1966; Stokes, 1972) and the 
Trignometric Collocation method (Samoilenko and Ronto, 1979). 

The primary objective of this work is to formulate the mathematical procedures for the 
analysis of periodic motion in nonlinear systems. The proposed procedure involves a coupling of 
the Trignometric Collocation method (TCM) with modal analysis techniques, thereby effecting a 
substantial reduction in the number of unknown quantities in the iterative part of the solution 
process. 
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MATHEMATICAL DEVELOPMENT 

The focus of this research is on the TCM that was developed and formalized by Ronto 
(Samoilenko and Ronto, 1979) with applicability to small order systems. Described below are the 
essential features of the procedure for both non-autonomous and autonomous systems. 

Trignometric Collocation Method 

Many engineering systems can be modelled by a set of n nonlinear ordinary differential 
equations of the non-autonomous form 


x = f ( x,t ) (1) 

where the RHS is continuous and periodic with a period T. It is required to determine a 
periodic solution x(t) of Eq. (1). It is assumed that the required solution can be approximated by a 
finite trigonometric series: 



where co is is the fundamental frequency. The unknown coefficients of the above series can be 
ordered into a vector, 

c i = (ao» a l» bj, a2, t>2, .... am, bm)T (3) 

corresponding to each variable xi. 

The collocation method essentially consists of substituting the assumed solution form, Eq. 
(2), into the system state, Eq. (1), and requiring that the equations be identically satisfied at a 
specified number of points, N. This gives rise to N*n nonlinear algebraic equations which must 

be solved to obtain the unknown coefficients. For a unique solution, the following inequality 
must be satisfied 


N£(2m+ 1) 


(4) 


Rigorous investigations of the applicability and foundation of the TCM have been carried out by 
Ronto (Samoilenko and Ronto, 1979), and only the formalism of the procedure is presented here. 


The state variables can be evaluated at the collocation points in terms of the unknown 
coefficients leading to the form: 


O • 





where, 


(5) 


o . 


x i 


( Xj (t o ), X. (tj), 



( 6 ) 


is a vector of values of the trignometric polynomial at the collocation points. The array T is an N * 
( 2m + 1 ) transformation matrix whose elements are defined as: 
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j = l 

j = 2, 4, ... 
j = 3, 5, ... 


(7) 


T.. 

y 


C: 


cos [ (i-1) jJt/N ] 
sin[(i-l)(j-l)Jt/N] 


The derivative of each state variable can be expressed in a trignometric series and is 
obtained by differentiation of Eq. (2). Hence, the following relation, which is similar to Eq. (5), is 
obtained: 


= col DCj (8) 

The array D is a (2m + 1 ) square tridiagonal matrix of the form: 


0 +1 

-1 0 

0 +2 

-2 0 


0 +m 
-m 0 


(9) 


and the elements are given by. 


D i,i + i = ' D i + i,i = i / 2 ’ i = 2 ’ 3 ’ 


., (2m) 


The requirement that the set of system state equations, Eq. (1), be satisfied exactly at N 
collocation points leads to N algebraic equations of the form: 


col D c. = fj(I c,^) (10) 

where fj is the vector of the ith function evaluated at the N collocation time points. Hence, the 
collocation process yields N*n nonlinear algebraic equations in the (2m + 1 ) * n unknown 
coefficients. These equations are then solved using a secant method from standard subroutine 
packages of IMSL. 

If the system state equations are autonomous , then Eq. (1) may be written as 


x = f(x) 


(ID 
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The frequency or frequencies of self-oscillation for such systems are unknown a priori. The 
analysis procedure is essentially the same as in the non-autonomous systems and leads to a set of 
nonlinear algebraic equations, 

©I D C j = f.(I c.) (12) 

The number of unknown quantities is increased by one, the unknown frequency, to [ (2m+l)*n + 
1 ], and the number of collocation points must satisfy the following inequality to assure a unique 
solution 


N 2 (2m + 1 ) + 1 (13) 

Hence, this situation is a case of non-linear least squares and cannot be solved by the secant 
method used for non-autonomous systems. An IMSL developed procedure, however, based on the 
Levinberg-Marquardt algorithm for nonliunear curve-fitting can be applied to the autonomous 
problem and has succesfully yielded satisfactory results, for several problems. 

Rotor System Equations 

The equations of motion for a typical multi-shaft flexible rotor system can be written in the 
second order form 


Mq + Sq + Kq = Q s + Q b (q b ,q b ) 

or equivalently in the first order form 
Ax + B x = F 


where 



and 



D 

m' 


-M 

o' 

— 



, B = 




M 

D 


0 

K 


(13) 

(14) 


(15) 


The linear forces of the system are included in the vector Q and the nonlinear component forces are 
included in the sparse vector Q b . A direct application of the TCM to a large order system such as 
Eq. (14) would almost always be computationally untenable. Thus, to obtain a mathematical model 
that is sufficiently small for the TCM to be effective, it is necessary to reduce the order of the 
original model. 


It is propsosed here to develop a procedure to analyze the periodic motion of large order 
structural systems with nonlinear supports or pseudo supports by using the TCM in conjunction 
with modal analysis. This algorithm will reduce the original problem to a set of nonlinear algebraic 
equations involving only the physical coordinates which are associated with the nonlinear 
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supports. This would normally result in a substantial reduction in order hopefully rendering the 
TCM computationally tenable. 

The connectivity between the nonlinear coordinates and the system displacement is 
specified by the connecticity matrix £ and can be written as 

q b (t) = S q(t) (16) 

where q b is the displacement vector associated with the coordinates of the nonlinear supports. The 
connectivity matrix S is a sparse matrix consisting mostly of zero elements and a few unity 
elements that ensure that the displacements in the nonlinear supports are identical with the 
displacements at the corresponding connection points of the linear system. Hence, Eq. (16) is a 
statement of geometric displacement compatibility for the system. 

Modal Analysis 

The A and E. arrays of Eq. (14) are not generally symmetric, thus both the right 
eigenvectors yj and adjoint left eigenvectors Zj must be evaluated for use in a modal expansion. 
These two sets of vectors satisfy the biorthoganality conditions 


Ay, - R, S. a) 

By, = -Wij b > 

where Rj is the system norm associated with the eigenvalue Xj. 


(17) 


With the state vector defined in Eqs. (15), the system eigenvectors are of the form 



where rj and lj are the right and left displacement eigenvectors associated with the physical 
coordinate vector q. The state response of Eq. (14) is represented by the modal expansion 


x = X^j < 19 ) 

i=l 

where Tjj is the ith modal coordinate. Substitution of Eq. (19) into Eq. (14) and premultiplication 
by ZjT , using the biorthogonality conditions of Eq. (17), yield the 2n equations 

li. - x.,11, = i-lT (Q S + Q b ) i = 1, 2 2n (20) 

1 

These equations are still coupled due to the nonlinear force vector Q b . For a large order system, it 
is not normally necessary nor is it feasible to retain all the modal information when determining the 
system steady state response. Usually only n lower modes are retained in the modal expansion of 
Eq. (19), thus there are correspondingly n equations in Eq. (20). 
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SOLUTION PROCEDURE 

Following the TCM procedure, periodic solution forms are assumed for the system 
physical coordinates, the nonlinear subset of the system physical coordinates, and the generalized 
coordinates in the modal expansion, i.e.. 


j=l 

a cos(oxt) + b sin(cor)J 

+ 

•V s 

II 

2 cos(cyt) + lx sin(co.t)| 

m 

"-H! + X 

j=l 

| 2 cos(co.t) + b\in(ci).t)j 


a) 


b) 


c) 


( 21 ) 


By choosing N equally spaced collocation points and evaluating the variables of Eq. (21) at these 
time points, the following set of relations is obtained 



= T Q 

a) 


Oqb 

= I £ b 

b) 

(22) 

°n 

= icn 

c) 



where the ith column of Q corresponds to the variable qj (t) evaluated at each of the N collocation 
points. The i* typical column of Q is defined by Eq. (3). Similar definitions apply for the arrays of 
Eqs. (22 b,c). 

The unknown coefficient arrays (£,£. b ,£n ) are dependent and are related through the 
geometric displacement relation of Eq. (16) and the modal expansion of Eq. (19). From Eqs. (22 
a,b) and (16), 


X £b = I £ £ t (23) 

and by utilizing the form of the system right vectors, Eq. (18), the modal expansion for the system 
physical coordinates may be written as 


q = X r i = r n (24) 

i=l 

Thus, from Eq. (22a) and (24) the following relation between physical coordinate and normal 
coordinate Fourier coefficients is obtained: 


I £ = I £H r T 


(25) 
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The substitution of this constraint relation into Eq. (23) gives 

£b = i T £ T ( 26 ) 

The next step in establishing the solution is to apply the TCM procedure to the set of 
pseudo modal equations, Eq. (20). Utilizing the equivalent form of Eqs. (5) and (8), Eq. (20) may 
be written as 


where: 


0) I E> cj^ - Xj X Cj 1 ! = °fi 


lf<Q s (t 1 ) + Q b (t,))| 


f = TT 

1 R. 


\ T 1 (Q(t N ) + Q b (t N )) 


(27) 


(28) 


The elements of the °f| vector are the RHS values of Eq. (20) evaluated at the collocation points, 
and are functions of the nonlinear displacements q b and velocities q b . 

Equation (27) can be rearranged to the form 

Cj 1 ! = [ (cofi-XiJ^a^I)- 1 ] I? °fi (29) 

i = 1, 2, .... n 

which is a typical column of the array of modal coordinate Fourier coefficients, £Tn. The 
combination of relations (29) with Eq. (26) results in a set of nonlinear algebraic equations in terms 
of Fourier coefficients for the physical coordinate subset q b .Thus, the size of the problem has 
been substantitally reduced and the location of a solution is computationally more feasible. 

The iterative procedure for estimating the Fourier coefficients £ b can be summarized in the 
following steps: 

1 . Choose a starting value for £ b . 

2. Compute °q b using Eq. (22 b). 

3. Compute Cl using Eqs. (19). 

4. Evaluate and update the value for C b using Eq. (26). 

5. Check convergence between steps 1. and 4. 

The procedure involves the solution of a set of nonlinear equations and its success depends upon 
the effectiveness of the numerical routine utilized. The optimization routine based on the secant 
method from the IMSL subroutine library proved to be very effective with convergence being 
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achieved for a wide variety of starting points. The error norm in the algebraic equations appears to 
be a reasonable measure of the solution accuracy. 

NUMERICAL EXAMPLES 

The results of analyzing three example systems are presented below. The first two are 
relatively small order systems and are analyzed directly in terms of physical coordinates. The third 
example is a larger order dual shaft rotor system that utilizes modal analysis in conjunction with the 
TCM. 

Journal - Hydrodynamic Bearing System 

Consider a single rotating journal on a hydrodynamic bearing as illustrated in Fig. 1. With 
reference to this figure, O is the bearing center, C is the geometric center of the journal, and c 
represents the radial bearing clearance. The mass center of the journal is assumed to be displaced 
from the geometric center by the eg offset e. During rotation this offset gives rise to a rotating 
unbalance force which is synchronous with journal spin frequency. The converging wedge that 
arises due to the eccentricity of the journal gives rise to a pressure field in the fluid film that 
supports the load. 

The nondimensionalized equations of motion of the journal assume the form: 

•• 

v = F cos<|> + F sin<|> + u cos ( t + (5 ) + g 

(30) 

w = F r sin<J> - F t cos<{> + u sin ( t + 13 ) 

In these equations, v and w represent the nondimensional displacement coordinates of the journal 
center with respect to a fixed reference frame. The quantities g and u represent the gravity and 
unbalance parameters, and F r and F t represent the radial and tangential fluid film force components 
acting on the journal. 

Using short bearing theory, Reynolds equation can be integrated to obtain closed form 
expressions for the plain journal bearing force components, e.g. (Holmes, 1960). Thus, 

F = p f ft ( 1 + 2 e 2 ) e + 2e 2 (l-2<fr) l 
r " L (l-e 2 ) 5/2 (1-e 2 ) 2 J 

(31) 


F = + B r — e e + 

1 L (1-e 2 ) 2 2 ( 1 - e 2 ) 5/2 “* 

where, e = ( v 2 + w 2 ) 1/2 , <(> = arctan(w/v) 

and B is a ^earing parameter that is dependent on the fluid viscosity, and geometry of the bearing. 

Clearly F r and F t are highly nonlinear functions of the response variables. Typically, the 
journal equations, Eqs. (30) are linearized about the static equilibrium position. The resulting linear 
response corresponds to an elliptical orbit centered at the equilibrium position. Application of the 
TCM to this problem can yield an orbit which is quite different the linearized response as displayed 
in Fig. 2. 
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An interesting fact revealed by the TCM is that the higher harmonics in the response of the 
journal may not be negligible as contrasted with many nonlinear problems. In the results presented, 
at least 8 harmonics were required to obtain close agreement between the TCM and numerical 
integration. In this case, many analytical procedures, such as the Describing Function Method, 
which neglect harmonics above the fundamental would not adequately describe the dynamics of the 
journal. 

Flow Induced Vibration 

Consider the problem of quasi-steady analysis of the transverse galloping of a long prism of 
square cross-section (Blevins, 1977; Parkinson and Smith, 1964).The equation of motion for the 
single degree-of-freedom oscillator is 


y" + 2 C y' + y = n u 2 c f (32) 

where U is the non-dimensional velocity of wind and y is the non-dimensional vibration 
displacement. The non-dimensional aerodynamic force coefficient Cf is obtained by experimental 

measurements in a wind tunnel and can be approximated by a polynomial in a, the angle of attack, 
or equivalently (y' / U). 


c 


f 


n 


AUy' - ^-(yf + 




(33) 


From a curve-fit to experimental values, A = 2.69, B = 168, C = 6270, D = 59,900 (for a 
Reynolds number = 22,300) ; n is a mass parameter, which is a function of the prism dimensions 

and the density of air, £ is the linear viscous damping coefficient 

The second-order nonlinear autonomous equation (32) has been shown to exhibit self- 
excited oscillations and an analysis by the method of averaging was carried out by (Parkinson, 
1964). It was found that the amplitude (Ai) vs wind velocity curves for various values of the 

damping coefficient collapse into a single curve if normalized by nAi / 2 £. 

The first harmonic amplitudes obtained by the application of TCM are shown in Fig. 3. 
As is evident from the figure, the response exhibits a hysteresis loop. A choice of different initial 
guesses helped the procedure converge to the multiple solution points. It is identical to the figure in 
Parkinson and Smith (1964) and indicates that the collocation procedure developed here is valid for 
problems with multiple solution points. 

Dual Shaft Rotor System 

The dual-shaft rotor system with configuration shown in Fig. 4 includes a nonlinear 
bearing at station 6 and excited by rotating unbalance in shaft 1 and a static side load at station 6. 
Rotor (1,2) is modelled as a (6,4) station, (24,16) degree - of - freedom, (5,3) element assembly 
with stations as indicated in Fig. 4. Detailed rotor configuration data is provided in (Nelson and 
Alam, 1983). The rotating assemblies are connected to a rigid foundation by linear bearings at 
stations 1 and 7 and are interconnected by a linear bearing between stations 4 and 10. A nonlinear 
bearing with cubic stiffness variation and linear viscous damping connects shaft 1 to the rigid 
foundation at station 6. The nonlinear bearing force components are given by the relations: 
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Fy - * ( kj r + 


( 34 ) 


. 3. v 

V )- - c v v 

F z - - ( k, r + kjr 3 )^. - c w w 

where r = ( v 2 + w 2 ) 1/2 , and kj = 50,000 lbf/in, k 3 = 50 *10 9 lbf / in 3 , and c v = c w = 20 lbf- 
s/in. 

The unbalance distribution of the rotating assembly consists of a single concentrated 
unbalance at station 2 with a eg eccentricity of 0.95 mils. In addition, a static side load acts on the 
system at station 6. Shaft 1 spins at 80,000 Rpm and shaft 2 co-rotates at 120,000 Rpm. 

It should be noted that the linear subsystem is not totally constrained. Thus, an "artificial support" 
is added at station 6 to eliminate a singularity. This influence is then subsequently removed from 
the model by subtracting its influence in the nonlinear forces. A value of 10,000 lbf / in at station 6 
was arbitrarily selected for this system. The nonlinear radial force versus displacement is shown in 
Fig. 5. The linear bearing stiffnesses are 150,000, 50,000, and 100,000 lbf / in at station 1, 4-10, 
and 7 respectively. 

Displacement orbits, as determined using the TCM procedure, for this system are plotted in Fig. 6 
for two stations and a side load of 100 lbf acting in the negative z direction. The orbit distortion 
clearly indicates the presence of higher harmonics in the response. 

CONCLUSIONS 

A numeric-analytic procedure based on the trignometric collocation method has been 
developed and implemented for estimating the periodic response of engineering systems. The 
procedure allows for estimating periodic forced response and for locating limit cycles of self- 
excited systems. A component mode synthesis strategy coupled with the TCM extends the method 
to large order system application. 

Three example analyses are presented. Two of small order in physical coordinates and the 
third on larger order using die modal strategy. Preliminary indications are that this method may be 
very effective in estimating the periodic response of both small and large order systems. Additional 
work is required to further test its generality, to handle systems with subharmonic response and, to 
ascertain the stability of the located periodic responses. Study on the speed and accuracy of the 
necessary computational work should also lead to improvement in the utility of the approach. 


1-23 



REFERENCES 


Atherton, D. F., 1982 , Nonlinear Control Engineering, Van Nostrand Reinhold: London. 

Blevins, R. D., 1977, "How Induced Vibration," Van Nostrand Reinhold, New York. 

Hagedom P., 1982, Nonlinear Oscillations, Clarendon Press: Oxford. 

Holmes, R., 1960, "Vibration of a Rigid Shaft on Short Sleeve Bearings," J. Mech. Engr. 
Sciences, Vol. 2, No. 4m pp. 337-341 

Nataraj,C., Nelson, H. D. and Arakere, N., 1985, "Effect of Coulomb Spline on Rotor Dynamic 
Response," Instability in Rotating Machinery, NASA CP 2409, pp. 225-233. 

Nayfeh, A. H., 1981, Introduction to Perturbation Techniques, John Wiley & Sons: New York. 

Nelson, H. D., and Alam, M., 1983, "Transient Response of Rotor Bearing Systems using 
Component Mode Synthesis: Part VI Blade Loss Response Spectrum," NASA Grant NAG3-6 
Report. 

Parkinson, G. V. and Smith, J. D., 1964, "The Square Prism as an Aeroelastic Nonlinear 
Oscillator," Quart. J. Mech. & App. Math., Vol. 17, pt. 2, pp. 225-239. 

Saito, S., 1985, "Calculation of Nonlinear Unbalance Response of Horizontal Jeffcott Rotors 
Supported by Ball Bearings with Radial Clearances," ASME J. of Vibration, Accoustics, Stress 
and Reliability in Design, Vol. 107, pp. 416-420. 

Samoilenko, A, M, and Ronto, N. I., 1979, Numerical- Analytical Methods of Investigating 
Periodic Solutions, Mir. Publishers: Moscow. 

Stokes, A., 1972, "On the Approximation of Nonlinear Oscillations," J. Differential Equations, 
Vol. 12, pp. 535-558. 

Urabe, M., 1965, "Galerkin's Procedure for Nonlinear Periodic Systems," Archives of Rational 
Analysis, Vol. 20, pp. 120-152. 

Urabe, M. and Reiter, N., 1966, "Numerical Computation of Nonlinear Forced Oscillations by 
Galerkin's Procedure," J. Mathematical Analysis Applications, Vol. 14, pp. 107-140. 

Yamauchi, S. 1983, "The Nonlinear Vibration of Flexible Rotors: First Report, Development of a 
New Analysis Technique," J. ASME, Vol. 49, No. 446, Series C, pp. 1862-1868. 


1-24 











